library(readr) # für den data import
library(dplyr) # für das datawrangling
data_star <-
readr::read_csv(
paste0(
"https://raw.githubusercontent.com/sammerk/",
"Zuerich_Statistics/master/data/star.csv"
)
)Statistische Modellierung mit R
Ein Workshop für die PH Zürich
Methodology: The Science before Statistics
Jede statistische Modellierung gewinnt an Aussagekraft, je umfassender sie die inhaltliche Fragestellung abzubilden im Stande ist. Um aus der riesigen Fülle an Optionen geeignet und zielgerichtet auswählen zu können sind die folgenden Unterscheidungen oft sehr hilfreich.
Erkenntnisinteressen
Ganz grundlegend kann a priori das Erkenntnisinteresse von Studien in die folgenden vier Kategorien unterschieden werden:
| Deskriptiv | Explorativ | Explanativ | Prädiktiv |
|---|---|---|---|
| populationsbeschreibend | hypothesengenerierend | hypothesenprüfend | Datenpunkte vorhersagend oder imputierend |
| Bei welchem Anteil 15-Jähriger in Deutschland handelt es sich um funktionale Analphabet:innen? | Was sind potentielle Ursachen für genderbezogene Disparitäten im Analphabetismus? | Sind 15-jährige Jungen häufiger Analphabeten als 15-jährige Mädchen? | Mit welchen Variablen können Schüler:innen at risk erfolgreich identifiziert werden? |
Gütekriterien wiss. Erkenntnis nach Campbell (1957)
Für ein erfolgreiches Studiendesign und die anschließende statistische Analyse ist es sehr wertvoll sich vorab über Schwerpunkte besonders gewünschter Aspekte wissenschaftlicher Güte Gedanken zu machen. Insbesondere über die Unterkriterien Methodischer Strenge:
- Konstruktvalidität (Inwiefern ist die Interpretation der Messwerte angemessen?)
- Interne Validität (Inwiefern sind Assoziationen von unabhängiger [beeinflussender] und abhängiger [beeinflusster] Variabler als kausale Effekte interpretierbar?)
- Externe Validität (Inwiefern können die Schlussfolgerungen der Studie verallgemeinert werden?)
- Statistische Validität (Wie robust und angemessen sind die verwendeten statistischen Verfahren?)
Deskriptiv- und Inferenzstatistik
Deskriptive Statistik beschreibt vorliegende Daten (z.B. mit Effektstärken), während Inferenzstatistik Aussagen über den die Daten generierenden Mechanismus trifft. Beide können »eher einfach« oder »hoch komplex« sein und oftmals stehen sie in einem synergetischen Verhältnis (siehe Abbildung 1).
Bayesianisches und Frequentistisches Schätzen und Testen
Eine sehr heuristische Klassifikation inferenzstatistischer Verfahren stellt die Unterscheidung von statistischer Schätzung und Testung dar:
(Inferenzstatistische) Schätzungen (estimation with quantified uncertainty) treffen anhand von Stichproben Aussagen über Parameter der Grundgesamtheit (Population) aus der die Stichprobe gezogen wurde.1
1 Bspw.: »Mit 96%iger Wahrscheinlichkeit liegt die Analphabetismusinzidienz von 15-Jährigen in Deutschland zwischen .08 und .12«
(Inferenzstatistische) Hypothesentests bewerten anhand von Stichprobendaten die Gültigkeit von Hypothesen in der Grundgesamtheit (Population) aus der die Stichprobe gezogen wurde.2
2 Bspw.: »Nimmt man an, dass sich die Analphabetismusinzidienz von 15-Jährigen in Deutschland zwischen 2021 und 2025 nicht geändert hat, beträgt die Wahrscheinlichkeit der vorliegenden Daten p = .032«
Diese beiden Verfahren können sowohl im Rahmen der frequentistischen Statistik als auch der bayesianischen Statistik angewendet werden. Die folgende Tabelle gibt einen Überblick über die wichtigsten Werkzeuge:
| Frequentistische Statistik | Bayesianische Statistik | |
|---|---|---|
| Parameterschätzung | Konfidenzintervalle | Posterior Distributions |
| Hypothesentest | p-Werte | Bayes Faktoren & ROPE-CrI Procedure |
Hypothesenarten
Bayesianische wie frequentistischen Hypothesentests können unterschiedliche Arten von Hypothesen zugrunde gelegt werden:
- Punkthypothesen setzen Parameter gleich einer reellen Zahl; etwa \(H_0\text{: } \delta = 0\)
- Äquivalenzhypothesen nehmen Parameter in einem reellen Intervall an; etwa \(H_0\text{: } \delta \not\in\ [-.3, .3]\)
- Informative Hypothesen nehmen eine Ordnungsrelation mehrerer Parameter an; etwa \(\mu_{\text{Baseline}} < \mu_{\text{Imaginary Pill}} < \mu_{\text{Blinded Placebo}}\) (Buergler u. a. 2023)
Die Art der (falsifizierten) Hypothese entscheidet wesentlich stärker über den Informationsgehalt eines Hypothesentests als die Entscheidung für das frequentistische oder bayesianische Paradigma (Hoijtink 2012).
Dies ist am leichtesten anhand der Nullhypothese nachvollziehbar. Wird etwa die Nullhypothese \(H_0\text{: } \delta = 0\) verworfen, wird entsprechend die Alternativhypothese \(H_A\text{: } \delta \neq 0\) angenommen. Diese enthält aber quasi keine Information, da sie nur mit einer einzigen Beobachtung (d = 0.000000 …) verworfen werden kann und im kritischen Rationalismus gilt, dass eine Aussage umso mehr Information enthält, umso leichter sie verworfen werden kann (Döring und Bortz 2016).
Äquivalenzhypothesen können sowohl frequentistisch (z.B. TOAST-Prozedur in R und JASP, Lakens 2017) wie bayesianisch (z.B. ROPE-Ansatz Kruschke 2015) getestet werden. Für das Testen informativer Hypothesen liegen bayesianische Methoden in (u.a.) JASP und R vor (z.B. {bain}, Gu u. a. 2019) sowie in den frequentistischen R-Paketen {restriktor} (Vanbrabant 2020) und {ic.infer} (Grömping 2010).
Grundlagen der Regressionsanalyse
Regressionsanalysen sind ein sehr mächtiges Werkzeug um Zusammenhänge oder Unterschiede zwischen/in Variablen zu modellieren. Innerhalb der Regressionsanalyse kann sowohl geschätzt als auch getestet werden und dies jeweils bayesianisch oder frequentistisch.
Die Grundidee der lineare Regression ist in diesem interaktiven Applet veranschaulicht. Eine Abhängige Variable \(y_i\) wird also in der Regressionsanlyse als normalverteilt mit bedingtem Erwartungswert \(b_0 + b_1*x_i\) und \(\sigma\) angenommen:
\[y_i \sim \mathcal{N}\left( b_0 + b_1*x_i, \sigma \right)\]
Datengrundlage
Das wollen wir in den berühmten Daten aus dem STAR-Experiment (Student/Teacher Achievement Ratio) illustrieren. In diesem Experiment wurden Schüler:innen und Lehrer:innen zufällig auf Klassen mit kleiner (13-17 Schüler:innen pro Lehrer:in) und großer (22-25 Schüler:innen pro Lehrer:in) Klassengröße zugeteilt. Zusätzlich gab es Klassen großer Größe, die von einer:m ausgebildeten Hilfslehrer:in unterstützt wurde. Die Schüler:innen wurden über mehrere Jahre getestet.
Wir können die Daten herunterladen oder direkt programmatisch importieren.
Der Satzsatz sieht wie folgt aus:
glimpse(data_star)Rows: 26,796
Columns: 18
$ id <dbl> 100017, 100028, 100045, 100045, 100045, 100064, 100070, 100070…
$ sch <dbl> 28, 52, 41, 41, 41, 64, 40, 40, 22, 22, 22, 17, 17, 3, 3, 3, 3…
$ gr <chr> "K", "K", "1", "2", "3", "K", "K", "1", "K", "1", "2", "K", "1…
$ cltype <chr> "small", "reg", "small", "small", "small", "small", "reg", "re…
$ hdeg <chr> "BS/BA", "MS/MA/MEd", "BS/BA", "MS/MA/MEd", "BS/BA", "BS/BA", …
$ clad <chr> "1", "1", "1", "APPR", "1", "1", "APPR", "1", "PROB", "1", "NO…
$ exp <dbl> 3, 12, 20, 15, 5, 19, 2, 5, 9, 9, 2, 20, 16, 8, 23, 21, 14, 6,…
$ trace <chr> "B", "W", "W", "B", "W", "W", "W", "W", "B", "B", "B", "W", "W…
$ read <dbl> 476, 410, 507, 575, 610, 430, 495, 629, 418, 524, 627, 421, 45…
$ math <dbl> 602, 444, 572, 620, 655, 484, 576, 592, 489, 567, 603, 459, 46…
$ ses <chr> "F", "N", "N", "N", "N", "N", "F", "F", "F", "F", "F", "N", "N…
$ schtype <chr> "inner", "suburb", "suburb", "suburb", "suburb", "rural", "rur…
$ sx <chr> "F", "F", "M", "M", "M", "F", "F", "F", "F", "F", "F", "M", "M…
$ eth <chr> "B", "W", "W", "W", "W", "W", "W", "W", "B", "B", "B", "W", "W…
$ birthq <chr> "1980:1", "1980:2", "1980:2", "1980:2", "1980:2", "1980:1", "1…
$ birthy <dbl> 1980, 1980, 1980, 1980, 1980, 1980, 1979, 1979, 1980, 1980, 19…
$ yrs <dbl> 0, 0, 1, 2, 3, 0, 0, 1, 0, 1, 2, 0, 1, 0, 2, 3, 0, 0, 0, 0, 1,…
$ tch <dbl> 478, 893, 698, 701, 706, 1102, 679, 684, 351, 359, 365, 267, 2…
Die Variablen dieses Datensatzes sind die folgenden:
id: a factor - student id numbersch: a factor - school id numbergr: grade - an ordered factor with levels K < 1 < 2 < 3cltype: class type - a factor with levels small, reg and reg+A. The last level indicates a regular class size with a teachers aide.hdeg: highest degree obtained by the teacher - an ordered factor with levels ASSOC < BS/BA < MS/MA/MEd < MA+ < Ed.S < Ed.D/Ph.Dclad: career ladder position of the teacher - a factor with levels NOT APPR PROB PEND 1 2 3exp: a numeric vector - the total number of years of experience of the teachertrace: teacher’s race - a factor with levels W, B, A, H, I and O representing white, black, Asian, Hispanic, Indian (Native American) and otherread: the student’s total reading scaled scoremath: the student’s total math scaled scoreses: socioeconomic status - a factor with levels F and N representing eligible for free lunches or not eligibleschtype: school type - a factor with levels inner, suburb, rural and urbansx: student’s sex - a factor with levels M Feth: student’s ethnicity - a factor with the same levels as tracebirthq: student’s birth quarter - an ordered factor with levels 1977:1 < … < 1982:2birthy: student’s birth year - an ordered factor with levels 1977:1982yrs: number of years of schooling for the student - a numeric version of the grade gr with Kindergarten represented as 0. This variable was generated from gr and does not allow for a student being retained.tch: a factor - teacher id number
Um nicht gleich mit Mehrebenenregression starten zu müssen, erstellen wir als erstes einen Subdatensatz, der nur eine:n Schüler:in pro Lehrer:in enthält
data_star_sub <-
data_star %>%
group_by(tch) %>%
sample_n(1) %>%
ungroup()»Effekte« der Klassenstufe
Einfache lineare Regression
Zunächst wollen wir »Effekte« der Klassenstufe modellieren. Ich empfehle dingend Daten immer erst zu visualisieren und dann statistisch zu modellieren.
library(ggplot2) # plotsWarning: package 'ggplot2' was built under R version 4.4.1
library(ggforce) # sina plots
ggplot(data_star_sub, aes(gr, math)) +
geom_jitter() +
theme_minimal()Warning: Removed 120 rows containing missing values or values outside the scale range
(`geom_point()`).
Die unabhängige Variable ist hier noch nicht intervallskaliert. Dies können wir wie folgt ändern.
# Das Schuljahr K zu 0 rekodieren
data_star_sub <-
data_star_sub %>%
mutate(
grade = case_when(
gr == "K" ~ 0,
gr == "1" ~ 1,
gr == "2" ~ 2,
gr == "3" ~ 3
)
)
# Anschließend nochmals visualisieren
ggplot(data_star_sub, aes(grade, math)) +
geom_jitter() +
theme_minimal()Warning: Removed 120 rows containing missing values or values outside the scale range
(`geom_point()`).
Dann kann eine erste Regressionsgerade modelliert werden.
Frequentistische Punktschätzung
Mit der Funktion lm() werden mithilfe kleinster Quadrate \(b_0\) und \(b_1\) geschätzt.
mod0 <- lm(math ~ grade, data = data_star_sub)
mod0
Call:
lm(formula = math ~ grade, data = data_star_sub)
Coefficients:
(Intercept) grade
485.16 46.64
Frequentistischer Nullhypothesensignifikanztest
Die Standardinferenzstatistik zu einem solchen Modell ist sind p-Werte für die Nullhypothesen \[H_0:\; b_0 = 0\] \[H_0:\; b_1 = 0\]
summary(mod0)
Call:
lm(formula = math ~ grade, data = data_star_sub)
Residuals:
Min 1Q Median 3Q Max
-197.162 -28.080 -1.162 25.920 142.559
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 485.162 1.980 245.04 <2e-16 ***
grade 46.639 1.083 43.08 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 42.55 on 1265 degrees of freedom
(120 observations deleted due to missingness)
Multiple R-squared: 0.5946, Adjusted R-squared: 0.5943
F-statistic: 1856 on 1 and 1265 DF, p-value: < 2.2e-16
Frequentistische Konfidenzintervalle
confint(mod0) 2.5 % 97.5 %
(Intercept) 481.27755 489.04613
grade 44.51531 48.76346
Bayesianische Posterior Distributions
library(MCMCpack) #für bayesianische Schätzung
mod1 <- MCMCregress(math ~ grade, data = data_star_sub)
plot(mod1)summary(mod1)
Iterations = 1001:11000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
(Intercept) 485.17 1.980 0.01980 0.01980
grade 46.64 1.095 0.01095 0.01095
sigma2 1813.72 72.286 0.72286 0.72286
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
(Intercept) 481.3 483.8 485.18 486.50 489.06
grade 44.5 45.9 46.62 47.36 48.78
sigma2 1677.0 1764.0 1811.87 1861.82 1959.01
Während die Densityplots den Posterior der Schätzung charakterisieren, erhält man mit summary() auch die 95%CrI-Intervalle.
Bayes Factor
Mit einem Bayes Factor testet man im Falle unseres Models wie wahrscheinlich die vorliegenden Daten sind geben ein Intercept-Only-Modell vs. unser Model mod1. Da die Funktion lmBF() im Paket {BayesFactor} allerdings Daten ohne Missings erwartet, müssen diese noch gefiltert werden.
library(BayesFactor)Warning: package 'Matrix' was built under R version 4.4.1
mod2 <- lmBF(
formula = math ~ grade,
data = data_star_sub %>%
filter(!is.na(math))
)Warning: data coerced from tibble to data frame
Standardisierte Regression
Die in mod0 erhaltenen Regressionskoeffizienten konnten wir nutzen um zu interpretieren welche Punktzahl durchschnittlich im Kindergarten vorliegt (\(b_0\)) und was der durchschnittliche Anstieg in einem Jahr ist (\(b_1\)). Diese Zahlen sind insbesondere dann gut als (nicht-standardisierte) Effektstärke interpretierbar, wenn man mit der Skala vertraut ist. Eine Alternative dazu ist die z-Standardisierung der abhängigen Variable. Um dies zu veranschaulichen betrachten wir Effekte des SES (free lunch) auf die Mathematikleistung in der dritten Klasse. Die Rohdataten sehen wie folgt aus:
data_star_sub %>%
filter(grade == 3 & !is.na(ses)) %>%
ggplot(aes(ses, math)) +
geom_violin() +
geom_sina() +
stat_summary(
fun.data = mean_sdl,
fun.args = list(mult = 1),
geom = "pointrange",
colour = "#1bbc9d"
) +
theme_minimal()Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_ydensity()`).
Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_sina()`).
Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_summary()`).
z-standardisieren wir die abhängige Variable bleiben die Verteilungsformen gleich, lediglich die Skalierung ändert sich.
data_star_sub %>%
filter(grade == 3 & !is.na(ses)) %>%
ggplot(aes(ses, scale(math))) +
geom_violin() +
geom_sina() +
stat_summary(
fun.data = mean_sdl,
fun.args = list(mult = 1),
geom = "pointrange",
colour = "#1bbc9d"
) +
theme_minimal()Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_ydensity()`).
Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_sina()`).
Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_summary()`).
Der folgende Code kodiert automatisch die ses Variable als 0, 1-Variable um.
mod3 <-
lm(
scale(math) ~ ses,
data = data_star_sub %>%
filter(grade == 3 & !is.na(ses))
)Die Koeffizienten passend sehr gut zur Abbildung:
\[ \operatorname{\widehat{scale(math)}} = -0.27 + 0.49(\operatorname{ses}_{\operatorname{N}}) \]
-0.27 ist der Mittelwert der Gruppe mit Free Lunch, 0.49 die Differenz der beiden Gruppenmittelwerte gemessen in Standardabweichungen und damit wie ein Cohen’s \(d\) interpretiertbar.
library(effectsize)Warning: package 'effectsize' was built under R version 4.4.1
cohens_d(
math ~ ses,
data = data_star_sub %>% filter(grade == 3 & !is.na(ses)),
pooled_sd = F
)Warning: Missing values detected. NAs dropped.
Cohen's d | 95% CI
--------------------------
-0.51 | [-0.74, -0.27]
- Estimated using un-pooled SD.
Multiple Regression
Die multiple Regression erweitert die Idee der einfachen lineare Regression auf mehrere unabhängige Variablen:
\[y_i \sim \mathcal{N}\left( b_0 + b_1*x_{1i} + b_2*x_{2i}, \sigma \right)\]
Informative Hypothesentests in der multiplen Regression
# Informative Hypothese
library(bain)
## ANOVA als Regression
data_star_sub <- data_star_sub %>%
mutate(classtype = ifelse(cltype == "reg+A",
"regA", cltype))
mod5 <- lm(math ~ classtype - 1, data = data_star_sub) # Intercept mit `-1` entfernt
mod5
Call:
lm(formula = math ~ classtype - 1, data = data_star_sub)
Coefficients:
classtypereg classtyperegA classtypesmall
544.3 555.4 558.3
bain(mod5, "classtypereg = classtyperegA = classtypesmall; classtypereg < classtyperegA < classtypesmall")Bayesian informative hypothesis testing for an object of class lm (continuous predictors):
Fit Com BF.u BF.c PMPa PMPb PMPc
H1 0.000 0.000 4.095 4.095 0.478 0.428 0.461
H2 0.734 0.164 4.477 14.064 0.522 0.468 0.504
Hu 0.104
Hc 0.266 0.836 0.318 0.036
Hypotheses:
H1: classtypereg=classtyperegA=classtypesmall
H2: classtypereg<classtyperegA<classtypesmall
Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement. PMPa contains the posterior model probabilities of the hypotheses specified. PMPb adds Hu, the unconstrained hypothesis. PMPc adds Hc, the complement of the union of the hypotheses specified.
Schätzung der eines Cohen’s d zwischen großer und kleiner Klasse in math
data_star %>%
filter(gr == 3 & !is.na(cltype) & cltype != "reg+A") %>%
cohens_d(
math ~ cltype,
data = .,
pooled_sd = F
)Warning: Missing values detected. NAs dropped.
Cohen's d | 95% CI
--------------------------
-0.17 | [-0.24, -0.11]
- Estimated using un-pooled SD.